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The ground state entropy of network source location problems is derived at both 
the replica symmetric level and one-step replica symmetry breaking level using the 
entropic cavity method. The recursive relations for the cavity probability and en- 
tropy are obtained in a probabilistic way by focusing on the change of the ground 
state size under the cavity iterations. The resulting entropic message passing inspired 
decimation and reinforcement algorithms identify the optimal location of sources in 
single instances of transportation networks. The traditional belief propagation with- 
out taking the entropic effect into account is also compared. We find that in the 
glassy phase the reinforcement algorithm is fastest while the entropic message pass- 
ing inspired decimation yields a value of the fraction of source nodes in the optimized 
network closest to the ground state value. 

Keywords: cavity and replica method, disordered systems (theory), message-passing algo- 
rithms, communication, supply and information networks 

I. INTRODUCTION 

Statistical physics methods and ideas inherited from studies of disorder systems [l|, 0] play 
an important role in providing theoretical insights and developing low complexity algorithms 
for combinatorial optimization problems or constraint satisfaction problems. These problems 
become subjects of interest across a variety of different disciplines such as computer science, 
discrete mathematics, statistical physics, engineering and computational biology jH- One 
archetype that involves both continuous and discrete variables is the network source location 
problem where we try to find optimal locations of sources in transportation networks with an 
optimal flow pattern. Studies on this source location problem [4t7| are practically relevant 
to network design and optimization [§1 and can have widespread applications in the field of 
operations research (9|. 

Transportation networks are composed of nodes with either a surplus or a deficiency of 
resources. How to distribute the resources and to choose some (deficient) nodes for resource 
provisions becomes an important problem in network optimization. In the source location 
problem, we are given a random distribution of resource providing nodes (surplus nodes) in 
a network and the objective is to optimize the actual location of additional source nodes 
so as to achieve the minimal installation cost of the new source providers as well as the 
corresponding transportation cost 0]. After the optimization, the remaining deficient nodes 
serve as consumer nodes and a network-wide satisfaction is achieved. Previous studies 
showed that the cavity energy functions with continuous variables can be decomposed into 
composite functions parameterized by their energy minima, such that the recursions of cavity 
energy functions can be converted to simple recursions of probabilities, which simplifies the 
analysis a lot 0,0,12, 10]. The resulting picture is a cascade of phase transitions in the glassy 



phase and the fraction of source nodes in the final optimized network shows discontinuous 
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transitions with increasing installation cost and this implies that different configurations 
of the source and consumer nodes will become energetically stable as the installation cost 
varies [H, 0|. Two of them are the singlet regime where the consumer nodes are isolated and 
the doublet regime where the consumer nodes can be paired or isolated. Even in the singlet 
regime, there exists glassy phase when the fraction of surplus nodes is below 0.25 for a 
random regular network with node degree equal to 3 Q. In the glassy phase, optimizing the 
locations of sources becomes algorithmically hard, which is most likely due to the emergence 
of long range frustration 11, 12]. Ref. 0] revealed that a random assignment of the bistable 
node with the same cavity energy in the consumer state or source state typically causes 
contradictions throughout the network. That is to say, the unfrozen nodes are strongly 
correlated in this region. Therefore, for those bistable nodes, we should consider their 
entropic effects and we expect that the entropic information extracted from the ground 
state can help to infer the correct assignment for each node. 

In this work, we assume that the cavity iterations 0], i.e., adding one node or one link 
between two nodes into the original system, will affect the size of the ground state (the 
number of optimal assignments or solutions), such that the recursive relations for the cavity 
probability and entropy can be derived in a probabilistic way. Hereafter, we identify an 
assignment to be optimal if it yields the minimal energy cost in a specific regime (e.g., 
the singlet or doublet connection pattern). This entropic cavity method has been used 
to compute the ground state entropy of the minimal vertex cover problem [l3[. Similar 
studies on the channel coding problem and budget-constrained auctions can be found in 



Refs. [14|, [15J. We emphasize here that the bistable node can either be in the source state 
or consumer state since both states share the same energy however, taking the entropic 
contribution into account, these two states may not have equal probabilities. Nonetheless, 
this entropic effect was ignored in previous studies 043]- Furthermore, by taking the zero 
temperature limit simply, the entropy term has vanishing contributions to the free energy 
and the ground state entropy can not be considered properly. In our current work, we 
shall demonstrate that not only can the entropic effect be used to evaluate the typical value 
of the ground state entropy as a function of the fraction of the surplus nodes (quenched 
disorder), but also can be used for algorithmic purposes on single random instances. We 
apply two kinds of decimation strategy to identify the optimal location of sources. One is 
the maximal decimation in which the most polarized nodes are fixed once the algorithm 
converges and then the network is simplified correspondingly. This strategy can be thought 
as a hard decimation since it is equivalent to adding an external field of infinite intensity to 
the decimated nodes EMI. The other is the reinforcement strategy which can be viewed 
as a sort of soft decimation. In this strategy, the cavity probability is strengthened or 
attenuated at each step by an external bias whose intensity is updated with a rate increasing 
with the run time (l9rm|. When the updating rate and the intensity of the external bias are 
correctly chosen, the external messages are able to drive the iterations towards some optimal 
assignment thanks to the fact that this soft decimation utilizes the global information of the 
network at each step. We find that in terms of computation complexity the reinforcement 
strategy is much more efficient than the hard decimation in the glassy phase. However, the 
entropic message passing inspired decimation is able to yield a value of the fraction of source 
nodes in the optimized network closest to the ground state value. 

The rest of this paper is organized as follows. The source location problem is defined in 
more detail in Sec. [Ill The analysis at both the replica symmetric (RS) level and one-step 
replica symmetry broken (1RSB) level by the entropic cavity method is presented in Sec. IIII1 



3 



In Sec. IIVI experimental studies of the proposed maximal decimation and reinforcement 
strategy are carried out. Conclusion and some future directions are given in Sec. [V] 



II. THE SOURCE LOCATION PROBLEM 

We consider a transportation network of N nodes in which resources are transported 
through the links, so that the resource demands of all nodes are satisfied. There are two 
kinds of nodes. The surplus nodes supply the resources and the deficient nodes have demands 
for resources. Respectively, <p s and 1 — <p s represent the fractions of surplus and deficient 
nodes in the network. The amount of resources supplied or demanded by a node i(i = 
1, . . . , N) is described by its capacity A* being positive and negative for surplus and deficient 
nodes respectively. To avoid considerations imposed by capacity limitations, we consider 
the generic case where the capacities of all deficient nodes are uniform and equal to — 1 
and those of all surplus nodes are uniform and equal to a sufficiently large real number 
A, so that resources from the surplus nodes are distributed to their neighbors with an 
optimal traffic pattern independent of the details of the capacity distribution. This scenario 
is typical in many modern applications of network traffic optimization, such as logistic 
networks and sensor networks, where the surplus nodes represent distribution centers in 
logistic networks or base stations in sensor networks. The locations of the surplus and 
deficient nodes are random, hence the capacity distribution of a node is given by p(A) = 
4> S S(A — A) + (1 — (j) s )5(A + 1). The network structure considered in this paper is a diluted 
network in which each node has the same degree C but the analysis below can be easily 
extended to the network with fluctuating degree (e.g., Erdos-Renyi random network). 

The source location problem is an important problem in dynamic network management. 
The task is to identify the locations of an appropriate number of deficient nodes and convert 
them to resource providers with capacity A ^> 1; these converted nodes are referred to as 
source nodes (the source nodes also include the surplus nodes when we consider the fraction 
of the source nodes in the final optimized network), whereas those nodes that remain as 
resource consumers are referred to as consumer nodes. The purpose is to reduce the overall 
transportation cost. On the other hand, the installation of source nodes has an installation 
cost, so that the total cost function is given by 

(ij) iev 

The first term is the total transportation cost, summed over the transportation cost of 
each link (ij) quadratic in the current which denotes the real-valued flow from 

node j to i. The quadratic function is chosen because of its convex property, which tends to 
balance the traffic load among the links; other convex functions can work equally well [IT 



241 ] . The flows are required to satisfy the non-negativity constraints of the final resources 
given by 

^ = Ai + ^^>0, (2) 

where £j is the final resource of node i and di denotes neighbors of node i. 

In the second term, u 2 /2 is the installation cost of a source node, where u has the same 
dimension as the flows. T> denotes the set of deficient nodes, and Sj = 1 if node i is assigned 
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to be a source node while for being a consumer node. Hence the source location problem is 
an optimization problem in the space of real-valued variables Xjj and the Boolean variables 

Si- 

It was found in Ref. |7| that the optimal choice of the source nodes can be obtained by 
assuming that all deficient nodes are consumer nodes, and then optimizing the nonlinear 
cost function in the space of the real-valued variables x^ given by [5] 

(ij) i€V 

where 0(x) = 1 if x > and Q(x) = otherwise. Then we identify those nodes with 
non-negative final resources as the consumer nodes, and assign those nodes with negative 
final resources to be the source nodes in the optimal solution of the source location problem. 

As the installation cost parameter u changes, different configurations of consumer and 
source nodes appear 0] • This leads to a cascade of phase transitions happen in the glassy 
phase with abrupt jumps of the fraction of source nodes in the optimized network. For 
instance, in the singlet regime where l/y/C < u < y/{C + 1)/[C(C — 1)], the consumer 
nodes are isolated. This is because in this range, the singly consuming state is always 
energetically more stable compared with its two neighboring phases such as all-source phase 
and doublet phase Q. The simple flow configuration in the singlet phase allows us to 
determine the optimal solution easily. The flow on a link into a consumer node is always 
1/C, and the flows on other links are Q. Hence the flow configuration is determined 
once the nodes are determined to be in the source state or the consumer state. This further 
simplifies the source location problem into a discrete-valued optimization problem. 

Based on the cavity method, a comparison between the energies of the cavity source 
state and the cavity consumer state led to a belief propagation algorithm 0, 0|. When 
the fraction of surplus nodes is sufficiently large, it converges satisfactorily and provides 
excellent agreement with the simulation results in terms of the fraction of source nodes Q. 
However, the algorithm was less satisfactory when the fraction of surplus nodes is not large. 
This is the regime where the cavity recursive relations become unstable towards fluctuations 
or, in the framework of the replica method, where the replica symmetric solution becomes 
unstable. 

One drawback of the previous analysis is that the degeneracy of the cavity source state 
and the cavity consumer state has been ignored. This was evident in Ref. [7j in which the 
cavity bistable state and the cavity source state were grouped together in deriving their 
recursion relations. Thus the entropy of the solution space was not considered properly. 
In this paper, we will compute the entropy of the ground state extending the previous 



works I7|,ll3[, restricting our attention to the singlet regime. We will also propose the entropic 
message passing inspired decimation algorithm and reinforcement strategy to identify the 
optimal location of sources in both the easy (replica symmetric) and hard (replica symmetry 
breaking) phases. 



III. THE ENTROPIC CAVITY METHOD 



In this section, we will present the entropic cavity method to compute the ground state 
entropy of the source location problem at the replica symmetric level and one-step replica 
symmetry breaking level. At the replica symmetric level, there exists a single ground state 
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and the clustering hypothesis that the correlation between any two randomly selected nodes 
is week becomes valid due to the sparse structure of the network. Following the work of 
Ref. [l3j, we focus on the change of the ground state size under the cavity iterations, and 
explicitly derive a closed set of equations involving two messages (cavity probability and 
cavity entropy. The replica symmetric ground state entropy can be evaluated from the fixed 
point of the recursive equations. Extending the analysis to one-step replica symmetry broken 
level is straightforward. When replica symmetry is broken, the single ground state would 
split into exponentially many states, which violates the clustering hypothesis. However, the 
weak correlation assumption is still satisfied in each state but the energy level crossings 
of states under the cavity iterations should be taken into account. All mean- field analysis 
presented here are restricted to the zero temperature limit which selects the ground state. 



Under the replica symmetric approximation, the joint state distribution of any two ran- 
domly chosen nodes from the large diluted network takes a factorized form making the 
derivation of a recursive relation feasible. Applying the cavity method to the source loca- 
tion problem, we consider the state of a node % in the absence of one of its neighbors j. In 
the singlet regime, it can be either in the cavity consumer state, or the cavity source state. 
In the cavity consumer state, the cavity energy of a node is lowest when a flow of 1/C enters 
it, whereas in the cavity source state, the cavity energy of the node is lowest when no flow 
enters it. 

Let us first define ipf^j as the cavity probability that node i is in a source state in the 
absence of node j. It can also be viewed as the message passing from the node % to node j. 
If V'f-s-j = 1, we say that node i is frozen to source state in the absence of node j. ipf^j = 
indicates node % should take the consumer state without node j. Otherwise, if ip'i^j G (0, 1), 
then the source and consumer states of node i are degenerate without node j. The value for 
the cavity probability V'f-^.j °f node i depends on the incoming cavity probabilities from its 
neighbors other than node j, which can be categorized into three cases 0. 



In the first depicted in figured] (a), all neighbors of node i have non-zero cavity 

probabilities {^j^}. In this case, the cavity state of node i must be a consumer. We 
assume that the number of optimal assignments before addition of node i is Qn-i where 
N is the number of nodes in the network. After the node addition in figure [T] (a), Qn-i 
should be reduced since node % is now in the consumer state and should remain as a singlet. 
Hence all its neighbors besides j should take source state and only those configurations with 
s k = l(k G di\j) in the ground state of the system with N — 1 nodes are valid after the 
addition of node i. di\j denotes the neighbors of node i except node j. In this case, V'iU? = 0) 
and the number of optimal assignments Qn = fijv-i ELeaAj rfk-H where the product comes 
from the week correlation assumption. The cavity entropy change is readily obtained as 



The second case (figured] (b)) gets a bit more involved. In this case, only one neighbor of 
node i, say node k, is frozen to the consumer state in the absence of node i, i.e., ipf.^ = 0. 



A. Replica symmetric analysis 
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FIG. 1: Entropic contribution due to node addition (adding node i and its adjacent links). 
Each node is connected to C = 3 nodes. The codes for the cavity probabilities are: shaded for 
"01 _>i G (0, 1] , black for ipf,^ = 0, white for ipf,^ = 1 and gray for ipf,^ G (0,1). The arrow 
shows the message passing direction, (a) None of the neighbors of node i is frozen to the consumer 
state in the absence of node i. (b) Only one of the neighbors of node i, say node k, is frozen to 
the consumer state without node i. After the addition of node i, it can take the source state or 
consumer state, (c) At least two of neighbors of node i are frozen to the consumer state in the 
absence of node i. 




(b) 



FIG. 2: Entropic contribution due to edge addition. The probabilistic meaning for each node is 
the same as that in figure [TJ (a) At most one end of the added link is frozen into consumer state 
before the link addition, (b) Both ends of the added link are frozen into the consumer state in the 
absence of the link. Then either node i or j should change its state after the link is added. 



According to figure [T] (b), the outcome is that the cavity source and consumer states of 
node i are degenerate. In Ref. 0], node i in this case was treated as the so-called bistable 
node. To consider the entropy change, we note that before the addition of node i, the 
number of optimal assignments in the ground state is Qn-i = Ylk'edk\i ^Pt'^k where 

is the number of optimal assignments in the ground state of the network without node 
k and i. Note that when node k was frozen to the consumer state, the same constraint 
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as the first case, namely, that ip%.'-+k > ^ or a ^ ^' ^ dk\i, was already imposed on node 
k. Following the discussion on the first case, this also implies that Qn-i can be simplified 
to £l N -i = £lN-2e ASk ~ >t where AS k _>i = J2k , €dk\i^°Si ! k'^k- The following two paragraphs 
analyze separately the possibilities that node i takes the source or consumer states. 

If node i takes source state as shown in the middle panel of figure Q] (b), node k needs 
not change its state, therefore the number of optimal assignments after node addition with 
Si — 1 is the same as flN-i- Hence we have fijvU=i = flN-2£ ASk ^ 1 ■ 

However, if node i chooses the consumer state, then the node k should change its state 
from consumer to source since we focus on the singlet regime where no paired consumer 
nodes are allowed. This change does not impose any further restrictions on the set of its 
neighbors dk\i, since the neighbors of a source node can either be sources or consumers. 
On the other hand, assigning node i to be in the consumer state restricts all its neighbors 
besides k to be in source states only. Consequently, the number of optimal assignments after 
node addition with Sj = is ^zvU=o = ^n-2 Yliedi\k j ^t-H where node i is fixed to consumer 
state (see the right panel of figure [1] (b)). 

Since node i has the choices to be in those ground state assignments that have Si = 
or Sj = 1, the total number of optimal assignments after addition of node i is Qn = 
f27v_2e A5fc ^ i + Qn-2 Yliedi\k j^t-H- The associated entropy change can be expressed as 

ASi^j = log -J^- = log 



1 + e~ As ^' Vf- 

ledi\k,j 



Qn-i 

At the same time, the cavity probability V'I-h is determined by 



(5) 





Si = l 




Si = l 




Si=0 



(6) 



The third case where at least two of incoming VI-h fo r node i vanish is presented in 
figured] (c). The added node i should take the source state Q, i.e., if)^ = 1- In figured] 
(c), both an d ?/f-n are equal to zero, thus the number of optimal assignments before 

the node addition is £l N -i = ^N-sUk'eakXi^k'^kHi'eaiXi^i'^i w here 1^-3 is the number 
of optimal assignments in the ground state of the network without node k,l and i. After 
node i is added, node i is frozen to source state and the neighbors k and / need not change 
their states, as a result, the number of optimal assignments Qn after addition of node i is 
identical to fijv-i- We conclude that the cavity entropy change for the third case is zero. 

Collecting these results, we can write the entropic message passing equations in a compact 
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form as 



^• =i - ( n v-iu + e t 1 - ^)]© n ^ 

,kedi\j J k€di\j \ledi\j,k 



l—ti 



X 



1 + e~ As ^ I 



ledi\k,j rl-H 



(7a) 



A^=e( n ^ e e t 1 - ^)] n ^ 

.k£di\j J k&di\j k£di\j \l€di\j,k 



X log 



1 + e ~ ASk ^ ip 

ledi\k,j 



s 

I— >i 



(7b) 



where a pair of messages (cavity probability and cavity entropy) are involved. In the above 
analysis, the added node is assumed to be deficient node. However, a finite fraction of surplus 
nodes with very large capacities are present in the transportation network as the quenched 
disorder. The addition of a surplus node is assumed to have no entropy contribution to 
the network and its full and cavity probabilities are always fixed to be 1 since it is frozen 
to the source state by definition. Adding a surplus node will yield different cavity energies 
depending on the states of its neighbors since the consumer neighbor will draw resources from 
its adjacent surplus node 0]. This leads to the recursive relations if)*^ = £a 15 a + ^L^U^-i 
and ASi^j = 8\ i ^AS^Ia^-i, where ip. 



* 

i->3 I 



_i and ASi^j\A i= -i are given by Eqs. (17a]) 
and ( l7bl) respectively. The entropy change AS; of adding node % can be computed using 
Eq. ( ITbl) where di\j is replaced by di as 

As, t = e ( n ^ ) e io s^ + e t 1 - e (^)i ( n ^ 

\kedi ) kedi kedi \l£di\k 



x log 



1 + e" 



n k 

l£di\k 



(8) 



To obtain the entropy density of the network, we need to add the link between two 
randomly selected nodes and consider the entropy change due to this link addition. This 
includes two C&SCS clS depicted in figure [2] (a) and (b) respectively. In the first case where at 
most one end of the link takes positive cavity probability, then after the link addition, those 
configurations where both node i and j take the consumer state should be excluded from the 
ground state whose size is denoted by Q, therefore, the number of the optimal assignments 
in the current ground state should be Q' = Q — Q(l — ^>?_^-)(l — Vj-w) with the entropy 
change AS^ = log [l — (1 — ^>f_^-)(l — ipj^)]. In the second case where both ends of the 
added link are frozen into the consumer state before the link addition, the number of optimal 
assignments in the ground state without the link is Q = J7jv-2 ILe9i\j -^i-h Tli'edjXi rft'-tj = 
QN-2e ASi ~^ j e ASj ~* i where is the number of optimal assignments in the ground state 

without node % and j. After the link addition, either node % or node j changes its state to 
the source state. The number of optimal assignments in the current ground state becomes 
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ft' = W\ s . =1 + 0'| s . =1 where ft'| a4=1 = Ili'eSjV = ^iV-2e A5j ^ and = 

'edi\j^i'^i — Si ^ j - Thus the entropy change due to the link addition in the 

second case is AS^ = log [e~ ASl ^ j + e 5j ' ->,i ]. To sum up, the entropy change due to the 
edge addition is written as 

AS m = 0(^ + ^) log [1 - (1 - - (g) 

+ [1 - Q^Uj + ^U)] log [e~ AS ^ + e" A5 ^] . 

The entropy density of source location problem in the transportation network can be 
evaluated in the Bethe approximation [3[ through 

s = < A ^> - | (A%)> (10) 

where (•) denotes both the disorder average and the average over the cavity message distribu- 
tion. We evaluate the entropy density by population dynamics algorithm A population 
of M pairs of (-0| AS^j) is used to approximate the joint distribution ■, AS^j) 

and its components are uniformly updated by the new computed ones according to Eq. (0). 
Usually, a number of iterations T — % are used to compute the entropy value with % it- 
erations for equilibration. Note that Eq. (ITUl) gives a self- averaging entropy value in the 
thermodynamic limit in the sense that the typical value of the entropy computed by the 
population dynamics algorithm should be consistent with that computed on a single large 
network j3|]. 

The source location problem has an interesting connection with the minimal vertex cover 



problem [13|, [25| if we assume the source nodes as the covered nodes, and the consumer nodes 
as the uncovered nodes. The difference is, we introduce the surplus (frozen covered) nodes 
as the quenched disorder in the source location problem where the network is a random 
regular network with fixed node degree and the total energy cost consists of the installation 
cost for all converted source nodes and the transportation cost induced by the consumer 
nodes. This specific energy cost makes the reweighting factor defined in Sec. IIIIBI totally 
different from that in the minimal vertex cover problem. 



B. One-step replica symmetry breaking analysis 

When replica symmetry is broken, the single ground state would break up into exponen- 
tially many ground states plus an even larger set of metastable states acting as the dynamical 
traps for any greedy search algorithm. In this case, one should take into account the reshuf- 
fling of free energies of different states when cavity iterations are performed, therefore, we 
write the replicated free energy $(?/) [26[ as 

e -yN3>{y) ^ 



Y,e- yNfa = j ' dfe N ™^ (11) 



where a indicates each state and £(/) is the complexity function counting states with given 
free energy density /. A saddle point analysis of Eq. (ITT]) gives $(y) = /* — S(/*)/y where 
/* is determined by y = dH(f)/df. The inverse pseudotemperature y allows us to weight 
differently the various states according to their free energy densities while the usual inverse 
temperature (3 = 1/T selects the energy of equilibrium configurations. Actually, Eq. (ITT]) 
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corresponds to a decomposition of the Gibbs measure [27|. Here, our analysis is restricted 
to the ground state (/3 tends to infinity), with / tends to the energy density e. For the 
current problem, the energetic complexity E(e) can be computed by the following Legendre 
transform 



y(e-$), 
dy ' 



(12a) 
(12b) 



where E(e) can be computed from the parametric plot of e(y) and £(?/) by varying the value 
of y. We remark here that the complexity function of the source location problem vanishes 
at finite y 0] and this was also observed in studies of the minimal vertex cover problem 25 . 
Therefore, the Parisi replica symmetry breaking parameter defined by m = y//3 vanishes in 
the zero temperature limit. The parameter m is used to select the size of the investigated 



ground state, as studied in Refs. [27-29[ to compute the entropic complexity curve for 
random constraint satisfaction problems where the ground state energy can be set to zero 
and the replicated free energy $ reaches the maximum at y — oo(in this case, m e [0, 1] 30] ). 

Due to the proliferation of pure states, we write the recursive equations for the joint 
distribution P^fy?^, AS^) at the 1RSB level as 

J kedi\j 



(13) 



where dP k ^i{il) s k ^ AS k ^i 

given by Eq. (!7al) and Eq. (17bl) respectively. The reweighting factor e _3/Ai/I ^ J takes into 
account the energy change due to the addition of node i and its adjacent edges except (ij). 
Note that m = makes the contribution of the entropy change in this term disappear. The 
cavity energy values have been computed in Ref. 0. It was found that the cavity energy 
AE^j is equal to u 2 /2 minus 7 for every cavity consumer state in its neighborhood set 
di\j, where 7 = u 2 /2 — 1/(2(7) which is non-negative in the singlet regime. Finally, is 
a normalization constant. To simplify the analysis, we parameterize the joint distribution 
according to the discussion in Sec. IIII Al as 



dr^idASk^Pk^t^ AS k . 



The function J 7 and J 7 are 



Pi 



where J d(p^ j (AS l ^. 



i-tjWi-tji AiSi-fj 



1 and pUi+pUi+P* 



with the RS case where the messages merely consist of the pair (ipf^j, AS^j 



(14) 

1. Compared 
for each 

directed edge, here the order parameter turns out to be a survey of these messages at the 
1RSB level; p}_>j, and p*_^j tell us the probability of picking up a pure state at random 
and finding that the cavity state of node % is source, consumer and free respectively. These 
three surveys enable us to handle the state-to-state fluctuations at the 1RSB level. With 
this parametric representation, a finite y survey propagation equation can be obtained for 
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the source location problem as 



Ai,-1 



n 



pI- 



PUi 



pU^+pI 



pI- 



k£di\j 

^2k£di\iPk^i eV1 Ylle8i\i,k \Pl- 



Pl 



Pk- 



(15a) 



(15b) 



pUj = S AuA + 8 Ai _x(l - pl+j - p*^). (15c) 



After the fixed point of Eq. (1151) is obtained, the replicated free energy can be computed via 

C 



(16) 



where the average is taken over the capacity distribution and the survey distribution, and 
the replicated free energy shift A$j due to node addition (and its C edges) and A$(y) due 
to link addition are, respectively, 



-yA$i = 5 A . _x log 



-yu 2 /2 



II K-* eW7 + Pl« + PU] + e^V 7 - 1) J] (pj^ + p*_ 



+ <W log ] [ [pL^ 7 + pU + p3U] 



-yA$ (ii) = log 



(pI+j + ?\* v )(/'.; - + Pi->i) + [i - (pUj + pUj)(:Pj -H + PjU*)] e?/7 



(17a) 
, (17b) 



where we have used the energy changes calculated in Tables 5 and 6 in Ref . 0] • The energetic 
complexity can then be computed using Eq. (TT2]) . In the glassy phase, E(e) typically increases 
from y = up to the maximal point forming the first non-physical convex part, yet, with 
further increase in y, decreases down to the zero point where the complexity vanishes at the 
ground state energy (y = y*). This second branch of the complexity curve is the physical 
concave part 0. Actually, the zero complexity corresponds to the maximum of the replicated 



free energy since E = y d y <& from the Legendre transform. To compute the ground state 
entropy at the 1RSB level, we should fix y = y*. 

To derive the formula for the ground state entropy at the 1RSB level, we first write the 
replicated free energy —y&(y, m) = S(e, s) — ye + ms keeping a finite value of m 28, 3(| and 



at the end of the derivation, we get the ground state entropy via s 



9{-y®(y,m)) 

dm 



\m=Q 



where 



y = y* determined by E(?/*) = 0. The 1RSB approximation of ground state entropy density 
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reads, 

/ [ASje-y**] \ C I [AS m e-y AE ^] 



\ [e-y^] j 2 \ [ e -y AE m] j ' (18a) 

[ASte-v**] = j J] d^(VU. AS fc ^)e-^ AS, ({^, AS fc ^}) , (18b) 

[e" yA ^] = S^e^ j J] M-^* 7 + vU + PU) + - 1) II + Pm) } 

I k&di kedi ) 



kedi 



(18c) 

[AS {ij) e-y AE w] = J dP^^AS^dP^^ (18d) 

[ e -^ W) ] = +pUj)(pU +pU) + e " 7 1 1 - (pUj+pUMpU+pU)} ■ 

(18e) 

Note that in the average [■], A Si and AS(y) are still given by RS equations (jHJ) and (j5J) 
respectively but the incoming messages to compute them should be sampled from the joint 
distribution Eq. (I14p and the reweighting factor is also included. In addition, (■) in the 
first term and second term in Eq. (I18al) denote the averages over the capacity distribution 
and survey distribution and this average can be easily done by the population dynamics 
algorithm [l| , which yields a a typical value of the entropy density. 

We evaluate the entropy numerically by population dynamics algorithm. At the 1RSB 
level, we require a population of N pairs of (p°_^,p*_^-) with an additional population of M. 
pairs of (ipf^j, AS^j) associated with each element of the first population. Another relevant 
quantity is the fraction of source nodes f s in the final optimized network. Its typical value can 
be evaluated as f s = ■h ^2 i=1 (ipf) at the RS level or 1RSB level. The average is taken over 
the disorder and the RS message or 1RSB survey distribution, and ipf is the full or marginal 
probability given by ip? = ^a^a + V^Ia^-i where Ui=-i is determined by Eq. (17a|) in which 
di\j is replaced by di. 

In the following studies, we choose u = 2/3 and C = 3 which fall within the range of 
the single regime 1 / \/~C < u < a/(C + 1)/[C(C — 1)] 0| • The entropy values of the ground 
state at both RS and 1RSB levels are given in figure [3J We use (T, 7o) = (5000, 1000) for 
RS computation and (T, 7o) = (1500,500) for 1RSB computation. Only the mean value 
of the entropy density is plotted since the error bar is much smaller than the symbol size. 
In the RS population dynamics, the extremely small value for certain cavity probabilities 
can be observed when is sufficiently low, and this is equivalent to the divergence of the 
evanescent cavity fields [l4, 3l[. However, if we take a proper cutoff e for the nearly vanishing 



probability, this situation can be circumvented [13j and the resulting entropy value seems 
to be insensitive to the cutoff value as long as e > 10~ 6 in the stable region of RS solution. 
For the 1RSB computation, we adopt e = 10~ 6 (smaller values of e do not change much the 
result). In figure [31 the entropy density first increases when <fi s decreases, then reaches a 
maximum followed by the decreasing trend as the fraction of surplus nodes (disorder) further 
decreases. At a large value of <f> s , most of nodes are surplus nodes, the entropy should take 
a small value; when (f> s gets close to one, most of nodes are deficient nodes, which makes 
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FIG. 3: (Color online) Ground state entropy density versus the fraction of surplus nodes <p s . The 
vertical dashed line indicates the dynamical transition below which the RS description becomes 
incorrect. The RS results are obtained using different cutoffs e = 10 -6 and 10 -8 . The inset gives the 
typical ground state fraction of source nodes in the final optimized network. The network parameter 
C = 3 and the population size for RS computation is TV = 10 5 , and (M,M) = (2048, 1024) for 
IRSB computation. 

the source locations much more constrained and yields a small value of entropy. Due to the 
concavity of the entropy function, there should exist some maximum point and this point 
can be identified by the mean field computation in figure [3J The maximum implies that at 
the corresponding <j) s ~ 0.40, there exist a largest number of optimal assignments satisfying 
the singlet connection pattern among all values of 4> s . 

At 4> s = 0.25, the replica symmetry solution starts to be unstable against IRSB pertur- 
bations [3]. This is confirmed by the fact that above this transition point, the complexity 
calculated by IRSB equations ffl~5]) does not depend on the value of y and is always zero 
while it has a finite value below the transition point. This transition point is also called 
the dynamical transition threshold below which the local search process is usually trapped 
by the metastable states whose entropy density is called the complexity. For the source 
location problem, the replicated free energy reaches its maximum at finite y and actually 
the maximum corresponds to the ground state free energy 0]. In the zero temperature limit, 
one can thus obtain the ground state energy. Within the current context, the ground state 
entropy value under the IRSB ansatz can also be computed improving the RS prediction. 
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We also show the fraction of source nodes in the ground state in the inset of figure [3] where 
the 1RSB ansatz predicts a higher value compared with the RS ansatz. We compute f s 
through the linear relationship f s — (e — 1/(2C) + 4> s u 2 /2)/^f since we have the identity 
(1 — f s )/{2C) + (f s — (p s )u 2 /2 = e. In this sense, the f s serves as a measure of the ground 
state energy density. These values will be compared with those obtained from simulations on 
single instances in Sec. IIVI As shown in the inset of figure El the f s decreases as <j) s decreases. 
Actually, one deficient node causes an installation cost [u 2 /2) when converted to a source 
node while it leads to an increment of 1/(2C) of the transportation cost when remaining 
as a consumer node. Notice that the consumer state is energetically favored since 7 > 
in the singlet regime. Consequently, on one hand, as <p s decreases, the number of deficient 
nodes increases, on the other hand, to yield the minimal energy cost, some deficient nodes 
need to remain as consumer nodes while consumer nodes are not allowed to be paired in the 
singlet regime. The competition of these two effects leads to the decreasing trend of f s with 
decreasing <p s . 



IV. NUMERICAL STUDIES ON SINGLE INSTANCES 

In this section, we apply the entropic message passing approach derived in Sec. Ill I Al to 
study single instances of networks. Two kinds of decimation strategy are used to identify 
the optimal source location. One strategy is the maximal decimation, and the other is the 
reinforcement strategy. The belief propagation inspired decimation [tJ is also compared. 
Hereafter, we use the short-hand notations BPD for belief propagation inspired decimation, 
EMPD for entropic message passing inspired decimation and EMPR for entropic message 
passing reinforcement algorithm. 



A. Maximal decimation strategy 

For the maximal decimation strategy, we fix a fraction fd of N t unfixed nodes with the 
largest full probability ipf to their most probable states once the algorithm converges at 
t-th sweep (fd = 2% in our simulations). Each sweep consists of a sequential update of the 
messages on all the edges of the network in a random order. When fdN t < 1, only one node 
with the largest bias is fixed. In the glassy phase, the algorithm typically fails to converge, 
and this is due to the building up of the long range correlations of different parts of the 
network. Therefore, we carry out the decimation according to the average full probability 
(over a number of sweeps) instead of the instantaneous value of the probability |32|. When 
a node, say i is fixed, the network is simplified by the following decimation procedure. If the 
node is fixed to the source state, it will send out the message Vf-H = 1 to all its neighbors 
j G di regardless of the later recursions. Otherwise, all neighbors of node i should take 
source state and are fixed at the same time since no paired consumer nodes are allowed in 
the singlet regime. Correspondingly, these fixed source nodes send out a constant message 
ipf_^.j = 1 to their neighbors. However, if at least one of neighbors of node i has been fixed to 
consumer in this case, a contradiction will be reported and we restart the algorithm. After 
the decimation, the recursion is carried out on the unfixed nodes of the network until N t = 
and an optimal assignment of source location is obtained. One can compute the marginal 

probability either based on the belief propagation Vf-+j — $Ai,A + ^Aj,-i 1 — ELeSAj ;"0fc-> 
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derived using the energetic cavity method Q, or according to the entropic message passing 
equations derived in Sec. IIII A[ The pseudocode of EMPD is given as follows: 

EMPD algorithm 

INPUT: the network with a fraction 4> s of surplus nodes; a maximal number of iterations T max ; a 
predefined precision 77. 

OUTPUT: one assignment in the singlet regime or 'probably no solutions'. 

1. Initialize randomly ^JL^-, ASi->j for all edges of the network except for the edges connected to 
the surplus nodes; 

2. for t = 1 to t = T max do: 

2.1 a sweep of random sequential update of the messages {tpf^j, AS^j} for all edges according 
to Eq. 0; 

2.2 if \tpUjit) - ipUjit —1)| < 77 and |AS , .^ i (t) - AS^i - 1)| < 77 for all edges and the 
iteration has converged: goto line 3; 

3. if t < Tmax, compute the full probability ijjf for unfixed nodes and decimate the network, else do 
the time average of the full probability over the later T max /2 sweeps and decimate the network; 

4. if an optimal assignment is found, return the solution and stop; else if no contradiction is found, 
continue the decimation procedure on the smaller network (goto line 2); else if a contradiction 
is meet, return 'probably no solutions' and stop. 



B. Soft decimation strategy 

The reinforcement strategy has been used to find solutions for random constraint sat- 



isfaction problems |19N21|. Here, we introduce an external bias fa for each node, and its 
intensity is updated with a probability increasing with the running time as 1 — t~ r . The 
external bias fa is updated as fa = it if ipf < 0.5 and 1 — it otherwise. The bias strength 
7r G [0, 0.5]. This implies that the current value of the cavity probability of one node will be 
enhanced by a factor if its full probability at the preceding iteration is biased towards the 
source state. Thus the only modification to the original recursive equation Eq. (jTJ) is 

& . = fig) 

^ 3 fa^Uj+(i-^)(i-r, 



where ^iUj is the reinforced cavity probability while V'iUj is the original one computed from 
Eq. (JTJ). Both 7r and r should be optimized so as to properly guide the iteration to converge 
to an optimal assignment. In the numerical simulation, we fix r = 0.1 although other choices 
are also able to find a solution, e.g., r = 0.3. The algorithm is precisely described as follows: 
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EMPR algorithm 

INPUT: the network with a fraction cf> s of surplus nodes; a maximal number of iterations T max ; two 
empirical parameters it and r. 

OUTPUT: one assignment in the singlet regime or 'probably no solutions'. 

1. Initialize randomly tpf^j , AS^j for all edges of the network except for the edges connected to 
the surplus nodes and initialize the external biases {/ij} at random; 

2. for t = 1 to t = T max do: 

2.1 a sweep of random sequential update of the messages {t/'fL^, A5j_^} for all edges according 
to Eq. (0) and Eq. <flbj> : 

2.2 update all the external biases with probability 1 — t~ r ; 

2.3 assign Sj = 1 if m > 0.5 and otherwise; if {s{} is an optimal assignment (no contradic- 
tion), return the solution and stop. 

3. if t = T rnax , return 'no solutions at current values of it and r'. 



The inference results using the above proposed algorithms on single instances of size 
N = 2000 are compared in figure HI As expected, the entropic message passing algorithm 
yields a lower value of f s than the belief propagation without taking the entropic effects 
into account. It should be mentioned that EMPR can find the optimal assignments faster 
than other strategies in the glassy regime, This is because, on one hand, the presence of an 
updating external bias manages to drive the evolution of the reinforced cavity probabilities 
towards the ground state where the satisfying assignment is maximally aligned with the 
externally imposed direction; on the other hand, the hard decimation fails to converge in 
this regime and the time-average of the marginal probabilities is required, which increases 
the time complexity of the algorithm. However, in the RS region, EMPR gets slow because 
the proper values of (n, r) could not be easily found, probably due to the presence of the 
many background messages sent out by the surplus nodes. On the other hand, the EMPD 
and BPD are typically convergent in this regime and thus fast to identify the optimal 
assignment. In the glassy regime, the estimated f s by message passing algorithms is higher 
than the ground state value predicted by the 1RSB theory. The physical interpretation 
is, in the glassy region, exponentially many metastable states act as dynamical arrests for 
various simple heuristics. Furthermore, the ground state predicted by the 1RSB theory 
can be unstable towards further steps of RSB or full RSB (an infinite hierarchy of nested 
states) jl, 0] , and the lower bound to the true dynamical threshold of f s is predicted by the 
Gardner value 33] above which the 1RSB metastable states become unstable jslj. However, 
the entropic message passing algorithm on single instances yields a lower value of f s , and 
particularly the EMPD can give a lowest value of f s among all message passing inspired 
algorithms in the glassy region. 



V. CONCLUSION 

We have considered the entropic effects in the ground state of the source location problem 
and derived the associated entropic message passing algorithms as an improvement over the 
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FIG. 4: (Color online) The fraction of source nodes in the optimized network (N = 2000) obtained 
by different algorithms. For BPD and EMPD, T max = 2000 and i] = 10~ 3 and additionally for 
EMPD, e = 10 -6 . T max = 10 4 for EMPR. The data point is the mean over 30 random regular 
networks with C = 3 and the error bar is also shown. The entropic message passing algorithms 
show lower values of f s close to the ground state values predicted by the 1RSB theory in the glassy 
regime. 



previous energetic computation of the problem. We assume that the cavity iterations will 
affect the size of the ground state and a self-consistent relation between a pair of messages 
(cavity probability and cavity entropy) can be obtained in a probabilistic way. Using the 
one-step replica symmetry breaking ansatz, the ground state entropy predicted by the RS 
solution is improved and yields a better approximation. Additionally, the predicted f s are 
checked by the simulations on single instances of transportation networks. Using the message 
passing inspired decimation algorithms, we find that the entropic information helps to lower 
the inference value of f s making it closer to the ground state value. 

As the installation cost parameter u increases, another important configuration of source 
and consumer nodes, namely the doublet regime will appear. Extension of the current 
method to this regime would be very interesting. In this case, to derive the recursive 
equations, we need to define two extra cavity messages, one for singly consuming state and 
the other for doubly consuming state. We leave this for future work. It is also of interest to 
extend the current analysis to the lattice glass models 35|, |36| where for example each site 
in the (finite dimensional or Bethe) lattice can have at most one particle and any particle 



18 



has at most a fixed number of occupied nearest neighbor. Other possible applications may 
be found in routing and path selection problems on sparse graphs 
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